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Abstract 

The separation of vectors by multigrid (MG) algorithms is applied to the study of 
convergence and to the prediction of the performance of MG algorithms. The separation 
operator for a two level cycle algorithm is derived. It is used to analyze the efficiency of the 
cycle when mixing of eigenvectors occurs. In particular cases the separation analysis reduces 
to Fourier type analysis. The separation operator of a two level cycle for a Schrodinger 
eigenvalue problem, is derived and analyzed in a Fourier basis. Separation analysis gives 
information on how to chose relaxations and inter-level transfers. Separation analysis is a 
tool for analyzing and designing algorithms, and for optimizing their performance. 
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1 Introduction 

The efficiency of multigrid (MG) techniques in solving large scale eigenvalue problems de- 
rived from discretizations of partial differential eigenvalue problems, was shown in multiple 
works, for example in [1] [3] [13] [16] [18] [20] [21] [27]. 

This work is motivated by the need of tools of analyzing, and designing robust and 
efficient MG eigenvalue solvers. These tools were needed for the algorithms presented in 
[4], [5] [6] [19] [24] and in the reports [7] [8] [9] [10]. The algorithms were applied to elec- 
tromagnetism and Schrodinger eigenvalue problems in 2D and 3D. Their efficiency resides 
in their eigenvector separation effectiveness, i.e., in how well they separate the subspace of 
desired eigenvectors from the remaining eigenvectors and how well they separate the eigen- 
vectors inside this subspace. These naturally hint to analyze how a procedure separates 
given eigenvectors, i.e., how it acts on the amplitudes of the eigenvectors, or more generally 
on the amplitudes of a basis of vectors. The action of an algorithm on the amplitudes of 
the vectors of a basis is called the separation operator. Separation factors, showing the 
amplification of one vector with respect to another one, can be used to define and analyze 
different convergence rates. In particular cases, Fourier analysis can be used to compute 
separation factors. Separation analysis does not reduce to convergence analysis only; for 
example, one may have divergence of relaxations but good separation, and separation de- 
termines the efficiency of the algorithm. Convergence results are obtained, which can be 
applied directly to show the convergence of single level or multilevel algorithms. The results 
are used for the analysis of a two level separation operator which defines the action of a 
two level algorithm on the amplitudes of given vectors. For eigenvalue problems, a major 
issue is the "mixing” of eigenvectors due to inter-level transfers. It is shown that better 
convergence is obtained if the mixing is small or when all eigenvectors which get mixed are 
treated together, e.g., in case of clusters. 

An example of separation analysis which reduces to a Fourier type analysis of the two 
level separation operator is presented. This is used in the design of MG eigenvector algo- 
rithms and in the prediction of their performance. For the Schrodinger eigenvalue problem 
with periodic boundary conditions, the analysis calculates accurately the convergence fac- 
tors. The matrix of the separation operator can be computed by a subroutine and the 
efficiency of the algorithm can be optimized analyzing the separation operator as a func- 
tion of several parameters such as relaxation types, iteration numbers. 

For MG techniques and more on MG eigenvalue algorithms we refer to [2] [3] [16]. An 
outline of general MG approaches related to separation and MG optimization is presented 
in [12] [25]. For eigenvalue algorithms and theory on algebraic eigenvalue problems we refer 
to [17] [23] [26] [14]. 
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2 The Separation Operator, Separation Factors and 
Convergence Rates 

A central goal is to analyze the action of an algorithm on the amplitudes of given vectors, 
for example on the eigenvectors of a given matrix. This action is defined by the separation 
operator. The separation operator can be used in the analysis of convergence and in the 
optimization of the algorithm. 

2.1 The Separation Operator and Separation Factors 

Assume that the n x n matrix U has independent columns Ui,...,U n . For example £/, may 
he Fourier components or the eigenvectors of a given matrix. Assume that an algorithm 
transforms the vector UB into the vector UB , where B and B denote vectors of dimension 
n. Thus B and B are the amplitudes of U\ , U n before and after the algorithm is applied. 
The algorithm defines a mapping 

WB = B (2.1) 

W is the mapping of the amplitudes of the vectors Ui under the action of the algorithm. 
The mapping W will be called the Separation Operator. The W may not be linear, e.g., 
for MG algorithms which employ projections. For particular MG cycles W is linear and 
will be computed in next sections. 

Consider further the case when W is an n x n diagonalizable matrix with normalized 
eigenvectors Ei,...,E n corresponding to the eigenvalues pi,...,p n . For simplicity assume 
that |/ii| > |/( 2 1 > ••• > |^n| > 0, although the further analysis can be performed in a 
similar way for the general case too. £j will be called the dominant eigenvector of W. 

For given vectors B = (6 X , .... b n ) T , B = WB = (6 X , ..., b n ) T define the separation 
factor of U t relative to Uj and to B by 

Si h B = (bi/h)/(bi/bj) (2.2) 

when this is meaningful, (i.e., when bj ^ 0, bj ^ 0, or define 0 and oo values when only one 
is 0). Implicitly, the depends on W also. 

It will be said that L r , is amplified by B if b t 0. Denote by e pm the p’th component of 
the eigenvector E m . Define the asymptotic separation factor of U % relative to Uj by: 

Sij = Ptlvi (2.3) 

where t = min{m : e, m ^ 0} and / = min{m : e.j m ^ 0}. The t and l in the definition of 
$ij, show the first (smallest index) eigenvectors E t and Ei which amplify U t . respectively 
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Uj. It will be also said that E t is the first eigenvector amplifying {/,. If E t is the eigenvector 
amplifying first both U, and Uj then 

sn = 1 (2.4) 

Sijjs, = 1 (2-5) 

2.2 Convergence of the Iterative Algorithm 

If an algorithm with separation operator W is iteratively applied starting with UB°, then 
the same results can be obtained by iterating W starting with B°. Thus the iteration of 
the algorithm reduces to a power iteration for W, just in the space of the amplitudes. The 
algorithm may be complicated but the power iteration is very simple (even for nonlinear 
W the same power iteration applies). The power iteration can be analyzed if W is known. 
Consider the representation of an algorithm by the iterations on W, in which the vector of 
amplitudes is normalized at each iteration (not necessary but to keep the norms bounded): 

Fourier Power Iterations 

B° = B 

For k — 1,2. ... do: 

1) B k = WB k ~ x 

2) B k = B k /\\B k \\ 

3) U k = UB k 

The following lemma shows towards what converge the Power Iterations and the algo- 
rithm. and shows the relative convergence rates. 

Lemma 1 In the Fourier Power Iterations algorithm, if F\ is the dominant eigenvector 
of W and if B° is not defective in E\ then: 

1) B k converges to Ej , 

2) U k converges to UE\. 

3) If B° is not defective in the first eigenvectors amplifying Ui and Uj then: 

s ij,B k s ij (2.6) 

Proof 1) is obvious since £j is the dominant eigenvector of W. Thus 2) holds. Let 
B° = £” =1 a m Em with di / 0. Then B k = J2m=i a mFmE m - Assume that E t and Ei 
are the first eigenvectors amplifying U t respectively Uj. The a t and a.{ are nonzero by the 
hypothesis 3). Then for sufficiently large k the denominator is not 0 and 

Sij.B* = (£™=f amj4 +1 e; m)(£m =( )/( (£" «m^ +1 e jm )(£” =< )) gt/ M = 
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In typical situations, e.g., in eigenvalue algorithms, it is desired that U k converges to 
l . Different convergence rates can be defined using W , its eigenvalues and eigenvectors, 
and the separation factors. Define the asymptotic convergence rate of U k towards U\ by 

Ci = sup oo max ll/sj^ ,B k I (2-7) 

£0 3=2.. .n 

where the B° in the sup is not defective in any of the eigenvectors. 

The following corollary is obtained directly: 

Corollary In the Fourier Power Iterations U k converges to U \ , for any B° nondefective 
in Ei. if and only if E\ = (1,0, ...,0) r . In case of convergence, the asymptotic convergence 
rate is ci = \p2\/\pi\- 

OBSERVATIONS 

The following two examples, when the convergence rate can be very good, motivate the 
above discussion: 

1) The inverse pow r er algorithm for eigenvalue problems: 

AU = C/A (2.8) 

where U are the eigenvectors of A associated to the eigenvalues of A = diag( Aj, ..., A n ). The 
inverse power algorithm iterates the operator (A — <J>/) _1 starting with U B° and normalizes 
the result at each step, for <l> a Aj . This corresponds to a Fourier Power Type algorithm 
where W = diag(l/(Xi — 4>)), (since {A — $I)~ l UB = UWB ~ UB), and with another 
normalization at step 2) (for which a similar result can be shown). The eigenvectors of W 
are Ei = U , and the eigenvalues are p, = 1 /( A ; — $). In this case the convergence rate is 
very good: |$ — Ai|/|4> — A 2 1 ~ 0; 

2) An MG Eigenvalue Cycle, where one may expect a similarly good convergence rate since 
the MG cycle can be viewed as an approximation to the inverse power iteration. Such an 
operator will be analyzed next, for a Two Level Cycle algorithm. 


3 Separation Analysis of the Two Level Cycle 
Algorithm 

This section derives and analyses the separation operator for a two level MG cycle algebraic 
algorithm. The cycle has the following useful features: 

1) The algorithm is algebraic and has no relation with any grid representation, neither the 
problems should be finer or coarser, the operators and transfers are general matrices; 
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2) The algorithm is simultaneous for several solutions, leading to a generalization of the 
separation setting from previous section; 

3) The simultaneous algorithm has advantages over the corresponding sequential algorithm, 
e.g.. it can incorporate on any level simultaneous separation techniques; 

4 ) The form of the algorithm is general so that the algorithm and its analysis can be used 
not only for eigenvalue problems but for system solvers too. 

In this cycle q solutions are treated simultaneously, thus the amplitudes matrix B has 
dimensions n x q. The names coarse level and fine level are used only to distinguish the 
two levels and because these names are used often in connection with MG algorithms. 

3.1 The Two Level Cycle Algorithm 

Suppose that A, U, A are n x n matrices, and that A\ U\ A' are m x m matrices such 
that: 


AU = U A (3.9) 

A'U' = U'lY (3.10) 

where U, U' are the eigenvectors of A, A' respectively, corresponding to the diagonal 
matrices of eigenvalues A and A'. In case of equal eigenvalues, the eigenvectors will be 
chosen to be linearly independent. Let P and R be two relaxation polynomials such that: 


P(A)U = UP{ A) (3.11) 

R{A')U' = U'R{iY) (3.12) 

There are given the transfer matrices J and 1 such that: 

JU = U'G (3.13) 

IU' = UF (3.14) 


where J and G are m x n matrices while / and F are n x m matrices. 

Assume given an n x q matrix, B , a diagonal q x q matrix $, and an n x q matrix 
U 1 = UB of initial solutions. 

Consider the following Two Level Cycle: 

Two Level Cycle (t/ 1 ,!/) 

Input V = UB 

1) Relax U 2 = P(A)V 

2) Compute the residual: 5 = U 2 $ — AU 2 
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3) Transfer the residual: S' = JS 

4) Transfer the solution: U n = JU 2 

5) Compute the FAS right hand side: T' = R(A')U n 4- JS 

6) Solve (or Relax) the coarse level equation: R{A')U' 2 = T' 

7) FAS - Correct the fine level solution: U 3 — U 2 + I(U' 2 — JU 2 ) 

8) Relax: 0 = P(A)U 3 
Output U 

OBSERVATIONS 

1) The step 6) can be considered as a relaxation: U 12 = R(A')~ 1 T' 

2) The algorithm can be continued on other levels in the same way. 

3 ) A projection can be introduced at certain steps, multiplying the solution by a q x q 
matrix E. Since E depends on the solutions, the algorithm would become nonlinear in B , 
making the next analysis more difficult. 

4) Different relaxations can be used at steps 1) and 8). 

5) The FAS (Full Approximation Scheme) transfers at steps 5) and 7) follow the scheme: 
the problem F'(U') = T' is an FAS transfer of the problem F(U) — T if T' = F'(JU) + 
J(T — F(U))\ and the corresponding FAS correction reads U — U + I{U' — JU), see [2]. 

3.2 The Two Level Cycle Separation Operator 

This subsection computes the separation operator W for the algebraic Two Level Cycle. 
The two level algorithm transforms the initial solution U l = U B into the final solution 
I = UB. Since the columns of U are independent, the B is uniquely determined for the 
given U , allowing to define as in previous section 
The Two Level Cycle Operator W by: 

WB = B (3.15) 

Denote by 7„ the n x n identity matrix. Assume that R(A') is invertible and denote its 
inverse by 7? _1 (.4'). The operator W can be computed directly from the above algorithm 
and from (3.11-3.14) and is given by: 

Theorem 

1 ) The Two Level Cycle Separation Operator W is defined by: 

WB = P( A)( P(A)B + FR~ 1 (A')G{P{i\)B^ - AP(A)B) ) (3.16) 

2) If q = 1 then 

W = P(A)( l n + FR-HA')G(^I n - A) )P( A) (3.17) 
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Proof Since U 1 = UB, the relations (3.11-3.14) imply: 


V 2 = P{A)U l = P(A)U B = UP{i\)B 
s = u 2 <f> - AU 2 = U(P(A)B* - AP(A)B) = UX 

where 

X = P(A)£$ - AP(A)f? 

The transfers give: 

U n = JU 2 = JUP(A)B = U'GP(A)B 
j,S = JUX = U'GX 

and solving the coarse level equation: 

R{A')U n = R{A')U n + JS 

U' 2 = U n + R' l {A!)JS = JU 2 + R~ 1 (A')JS = JU 2 + U'R~ l {A')GX 
Then the FAS - Correction reads: 

U 3 = U 2 + I{U' 2 - JU 2 ) - U 2 + IU'R~ l (i\')GX 
U 3 = U(P{A)B + FR~ 1 (A')GX) 

and the final relaxation implies: 

U = P{A)U 3 = UP{A)(P{A)B + FR-\A')GX) = UB 

Hence: 

B = P{A){P(A)B + FR-\A')GX) 

Substituting A’, the point 1) of the theorem is obtained: 

WB = B = P{ A)( P(A)B + FR- l {A')G(P{A)B<!> - AP{A)B) ) 

For the case q = 1, $ is a scalar and commutes with P(A)B thus point 2) results. 
Iterating the Two Level Cycle the following algorithm is obtained: 

Two Level Iterations 
Input U 1 = UB , $ 

for k = 1 , 2 ,... do: 

1) Two Level Cycle (U k ,U k+1 ) 


(3.18) 

(3.19) 

(3.20) 

(3.21) 

(3.22) 

(3.23) 

(3.24) 

(3.25) 

(3.26) 

( 3 . 27 ) 

(3.28) 

(3.29) 



2) Normalize the columns of U k + l 


OBSERVATIONS 

1) For the Two Level Iterations algorithm holds a Fourier Power Iterations algorithm as 
described in section (2). 

2 ) As in the Corollary of Lemma 1 , for q = L the convergence criteria is obtained: 
Lemma 2 In the Two Level Iterations , U k converges to f/ 1? for any initial U l nondefec- 
tive in U i, if and only ifW has the dominant eigenvector (1,0,.. .,0). The single difference 
consists in the normalization of B k by B k — B k /c k where c k is a constant for normalizing 
the solutions. 

3) Denoting B = (&i, b n ) T , the k component of the vector WB is: 

n m 

(WB) k = P 2 (\ k )b k + P( A fc ) X] Z(Fk,R- l (K)G ij )(^ - Aj)P{\j)bj (3.30) 

j=l i = 1 

This formula will be used in the following sections. 


3.3 Separation Factors in case of Mixing 

In this paragraph the separation factors are analyzed in a case when the first two vectors 
of U are mixed by IT during the transfers, but are not mixed with other vectors, i.e., the F 
and G have the structures from (3.31, 3.32). It will be shown that the mixing damages the 
separation, nevertheless the algorithm may be efficient for proper choices of relaxations. 
Assume that q = 1 , B = ( 6 j, 62 , 0 , ...O ) 7 and: 


G = 


/ 1 
G 21 
0 


G 12 
1 
0 


Vo 0 


0...0 \ 
0...0 


F 


( 1 
Fn 
0 

V 0 


F 12 

1 

0 


0...0 \ 
0...0 


(3.31) 


(3.32) 
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The first two components of W B will be: 


(W B)\ = P 2 ( Ai)6i + P(Aj)((/? X (A [) + F 12 R 1 (A' 2 )G2i)($ — Ai)P(Ai)6i + 

+(i?- 1 (A / 1 )G’ 12 + FuR-'i A'))($ - A 2 )P(A 2 )6 2 (3.33) 

(WBh = F 2 (A 2 )6 2 + P(A 2 )((F 21 F“ 1 ( a;) + p- 1 (A')G 21 )($ - A 1 )P(A 1 )6 1 + 

+(F 2 ijR _, (Aj)G 12 + R~H A' ))($ - A 2 )P(A 2 )6 2 (3.34) 


Assume that the Richardson relaxation is used for A: 

P(A) = I + u{A - $) 


(3.35) 


such that 

F(Afc) = 1 + a;(A/t — 4>) ~ 1 

and the inverse power iteration is used for A'\ 


R(A') = A' -<t>, R{A k ) = A' k -$ 


(3.36) 

(3.37) 


To simplify the analysis it is assumed that: 

$ « A, (3.38) 

P{\i) « 1 ( 3 . 39 ) 

($-A 2 )/(A' 2 -$)*-l ( 3 . 40 ) 

A) ss Aj (3.41) 

The next approximation is obtained for the first two components (the rest of the matrix is 
not relevant and ignored further): 


Where 


WB « 


1 7 - P\2 
0 7F21 



7 — £12(4* — a 2 )/(a; — $) 


( 3 . 42 ) 


(3.43) 


If G 12 is not zero then 7 may be large since A[ — 4 > is small. The two eigenvalues of 
the above W are pi = 1 and p 2 = 7F21 corresponding to the eigenvectors E\ = (l,0) r 
and E 2 = ((—7 + Fi 2 )/(1 — 7F21), l) r . If 1 > |/t 2 | then iterating the two level cycle, the 
solutions UB k wdll converge to U\, by Lemma 2 . If 1 < |p 2 | then the iterations wdll converge 
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to (Fi,F 2 )F 2 . In the case 1 < |p 2 |, the separation can be improved by: 

1) improving the relaxation or the transfers (e.g., preconditioning the transfers or using 
higher order transfers), these leading to 1 > |p 2 |; 

2) treating simultaneously the vectors Ui and (/ 2 and separating them by a Rayleigh-Ritz 
type projection as shown for example in [5] [8] and analyzed in [10]. 

The asymptotic convergence rate in the case 1 > |// 2 | is 

|^2/pi| = |p 2 | < 1 (3.44) 

which can be good, e.g., if Gi 2 F 2 i is small, but can be close to 1 if 7 is large. Nevertheless 
the convergence in the first few cycles may be very good as shown by the relative separation 
factor. If at the beginning &i = 6 2 = 1 and 7 is large, then s 12 ,b = (1 + 7 — F 12 )/(7F 2 i) ss 
I/F 21 . In this case the separation factor of U\ relative to U 2 is approximately 1/F 21 which 
may be very good. A better separation is obtained if Gu ~ 0, when 7 « 0 thus the 
separation factor is (1 — Fi 2 )/(7Fi 2 ) « 00. In the latter case, when low mixing appears 
due to G'i 2 . the cycle separates well the first eigenvector from the second one in the first 
iteration. An important observation is that in this case the role of the relaxation by P is 
negligible, the good separation being due mainly to the good transfer G\ 2 ~ 0 and to the 
good approximation of the eigenvalue Ai ~ 

3.4 Separation Factors for Decoupled Components 

To show the good efficiency of the MG cycle in the case when the eigenvectors are not 
mixed during the transfers, the previous example in which the mixing coefficients are 0 is 
considered: 

G 12 — G' 2 1 = F\2 = F 2 i = 0 (3.45) 

A more accurate estimate of the relative separation factor for bi — & 2 = 1 and by (3.33,3.34) 
is: 

[1 + <* - A,)/<a; - #)]/[l + (* - A 2 )/(A' - *)] = (3.46) 

which is large in the following assumptions which are met especially in MG algorithms 
where 4> is obtained from coarser levels by an FMG algorithm [2] [5] [6]: 

A 2 - A 2 « 0 (3.47) 

(a; - ao/ca; - « 1 

|A' - $| > 0 
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(3.48) 

(3.49) 



OBSERVATION 

The efficiency of the two level cycle may be very good also in case of close eigenvalues, 
when $ approximates well the eigenvalue Ai and when the two levels have close eigenvalues 
A 2 ~ A 2 . This is an explanation of the highly accurate separation obtained in numerical 
tests for problems presenting very close eigenvalues [4] [5] [6] [8]. The good separation 
obtained by the MG cycle may also explain why projections were not required on fine 
levels even for clusters of very close eigenvalues [4] [5] [6] [8]. 


4 Example of Fourier Analysis of the Two Level 
Separation Operator 

In some cases, especially for discretizations of partial differential equations on regular grids, 
the two level separation operator can be computed and analyzed using for U Fourier compo- 
nents. as shown in the next example. This provides insight for the design of the algorithms 
and to the prediction of the algorithms performance. 


4.1 The Two Level Separation Operator 


The two level separation operator IT for the algebraic two level cycle, in the common 
assumption n = 2m, is presented next. For n = 2 m the matrices F and G have the forms 
F = (Fi,F 2 ) T , G = (Gq , (* 2 ). with m x m matrices F,, F 2 , G,, G 2 . Further, Id will 
denote the identity matrix. It is assumed q = 1, so W has the form (3.17). Denote 


A = 



0 

A 2 


where A,, A 2 are m x m diagonal matrices. In this case (3.17) provides: 

w = ( H n IT, 2 
{ W 21 W 22 

with 


lT n = P(A 1 )(/d+F 1 i?" 1 (A , )G' 1 ($/d-A 1 ))P(A 1 ) 
IT 22 = P(A 2 )(/d + F 2 R~ l ( A r )G 2 ($Id-A 2 ))P( A 2 ) 
II j 2 = P{Ai)FiR~ l (A')G 2 ($Id — A 2 ))P(A 2 ) 

14 21 = P(\ 2 )F 2 fF i (\')G\(d>Id - A, ))F(A,) 


(4.50) 


(4.51) 


(4.52) 

(4.53) 

(4.54) 

(4.55) 
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4.2 The Two Level Separation Operator for the Schrodinger 
Eigenvalue Problem 

The two level separation operator for the Schrodinger Eigenvalue Problem in 1 -D, with 
periodic boundary conditions, is derived further. For 2 -D and 3-D similar derivations hold. 
The problem is: 

AU = U A (4.56) 

on the interval [0, 27 t] . The operators A and A! are the discretizations of the Laplacian, 

with the stencil ( 1—2 1 )/h 2 , on two grids with 2 m (respective m) equally spaced points. 

The complex eigenvectors and the corresponding eigenvalues are considered: 

Uk{l) = exp(i2rkl/2m). X k = ~ — (l — cos{2irk/2m)), k = 0 2m — 1 (4.57) 

/>“ 

U'kti) = exp{i2irkj/m), X' k = ~ cos(2tt k/m)), k = 0,...,m - 1 (4.58) 

The eigenvalues are negative, ranging from A 0 = 0 to the size of \ m = —4 /h 2 . where 
h = 2ir/(2m) is the fine mesh size. The eigenvalues satisfy 

Xm — k Xfc — A 2771 — k 0 ,..., TTl 1 (4.59) 

The real and imaginary parts of the complex eigenvectors are the real eigenvectors: 

U r k (l) = cos(i2vkl/2m), k = 0, ...,2m — 1 (4.60) 

U k (l) — sin(i2wkl /2m), k = l,...,2m — 1, k ^ m (4.61) 

where i ' T 2rn _ k — and U' 2 m _ k = —U' m . The analysis can be done for the real eigenvectors 
but it- is more convenient computationally to use the complex ones. The real eigenvectors 
can be used in programs which build W and analyse the separation. 

The transfer J is the full weighting operator with stencil (12 1 )/4. / is the linear 
interpolation operator. Then the matrices F\, F 2 , G\. G 2 are diagonal and all elements of 


F and G are zero except for k = 0, .... m — 1: 

G kk = (1 + cos(7rfc/m))/2 (4.62) 

G k , k+m = (1 - cos(7rfc/m))/2 (4.63) 

F kk = (1 + cos(7rA:/m))/2 (4.64) 

F k+m ^ k = (1 - cos(nk/m))/2 (4.65) 

since the transfers imply for k = 0, ..., m — 1, j = 0, ..., m — 1: 

(f k(2j — 1) + 2V k (2 j) + V k (2j + l))/4 = V k (j)( 1 + cos(?rA’/m))/2 (4.66) 
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{U k+m (2j - 1) + 2U k+ m(2j) + u k+m {2j + l))/4 = U' k (j){l - cos(jrfc/m))/2 (4.67) 


Uk(2j) = U' k (j) (4-68) 

U k (2j + 1) = (U' k (j) + U' k (j + l))/2 (4.69) 

thus 

JUk = U k G kk , A: = 0, .... m — 1 (4.70) 

JU k + m = U' k G k , k+m , fc = 0, 1 (4-71) 

IU' k — F kk L' k 4- F k +m. k [ J k+m 1 k = 0, ...,m — 1 (4.72) 

Then the elements of the operator IT become: 

ITn = dfa<7(.P(A fc )(l + Ffcfci? -1 (A'jGfcfc(<I> - A fc ))P(A fc )) (4.73) 

VT 22 = diag(P(\ k+m ){l + F k+m , k R~ 1 (\' k )G k , k+m {$ - \ k + m ))P{\ k + m )) (4.74) 

IT 12 = dtag(P(X k )F kk R-^X' k )G kMm (^ - X k+m )P(X k+m )) (4.75) 

IT 2 i = diag(P(\ k+n )F k+m , k ir'(\' k )G k , k {* ~ A fc )P(A fc )) (4.76) 


4.3 Analysis and Optimization of the Two Level Cycle 

The TT’ operator (4.73-4.76) is used to analyse and optimize the two level cycle algorithm. 
R is taken the shifted inverse power operator 


£0 

1 

II 

1 

3 

(4.77) 

The relaxation operator is taken 


P{A) = Id + w(A-$) 

(4.78) 

then 


P(X k ) = 1 + u?{\ k — $) 

(4.79) 

Assume that 


U.* — r^/jA maa .| 0 /? /4 

(4.80) 

where a will be chosen in a convenient way. To analyze the relaxation 
asymptotic separation factor of relative to Uj and to the relaxation is 

note that the 

_ |1 + u;(Ai — A)| _ ||A max | + q(A 1 — A) 1 _ 
|1 + u;(Aj — A)| ||A mar | + ©(Aj — A ) | 

(4.81) 

a( Aj — Aj ) . 

|^mar| + a ( ~ 

(4.82) 
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It follows that if A, ss A_, then s tJ ~ 1 thus the relaxation will be very slow in separating 
U, from Uj. If A, is close to Aj and A j is close to \ max then s tJ & |1 — a/ (1 + a)|. Thus the 
relaxation is more efficient in separating far clusters. Within clusters another separation 
technique, like a Generalized Rayleigh Ritz projection, or MG projections, or an MG cycle is 
required, see for example [5] [8] [10]. An essential property of the relaxation is the damping 
of errors coming from eigenvectors which are not well represented on coarser levels. This 
is called the smoothing property of the relaxation. If U, is a smooth vector and Uj is a 
component of the error, then Uj is damped by each relaxation with the separation factor 
s tJ . 

Next computation finds an a such that the separation factors s 0 h > > 1 for all 

nonsmooth components Uk , assumed A = A 0 = 0. Recall from (4.57) that the nonsmooth 


components Uk have the frequencies tt/2 < 6 = tt kjm < 3ir/2 corresponding to m/2 < 
k < 3m/2 and to |A max | = |A m | = 4 /h 2 . Then 

l-Sojtl = |2/(2-a(l - co.s(7rA-/m))| (4.83) 

The extreme values of |so*-| are obtained for 2irk/2m — 9 = 7r/ 2 and 7r, thus k = m/2 
respective m, for which |s 0m / 2 | = (2/(2 — a)| and |s 0m | = 1 1 /( 1 — a)|. Both 7r/2 and k 
components can be separated by relaxation from Uo with the factor |.s 0 to / 2 | = |som| = 3 if 

q = 4/3 (4.84) 


Moreover, one ralaxation will damp all high frequency components by a factor larger than 
1/3. If a multilevel cycle damps well the smooth components of the errors, on coarse lev- 
els where these components can be well represented, and does not amplify the oscillatory 
components, then the cycle should have a separation factor of (1/3) M , for p relaxations per- 
formed on the fine level. This factor may not be obtained in case of mixing of eigenvectors, 
and in the case when not all smooth components are damped well, e.g., the frequencies 
close to the frequencies of the desired eigenvectors, are damped slowly by the relaxation. 


OBSERVATIONS 

1) The WB is 


WB = W 


B i 

B 2 


WnBi W u B 2 \ = f Br \ = - 

W' 2 \Bi W22B2 J y B2 J 


(4.85) 


Thus if U 21 is not small enough then the large coefficients of the smooth components, Bj , 
of order 1. will imply large coefficients B 2 of the nonsmooth errors. The factors of W 2 i 
which can be improved in (4.76) are 


P (Afc+TTj ) Fk+m.k 


(4.86) 
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This can be done by changing the value of a, say close to 2, to damp more efficiently the 
frequencies close to x/2, 3 7T / 2 . which are multiplied with 1/2 by transfers. 

2 ) The relaxation can be improved also using a combination of relaxations with different 
values of a which will damp different frequencies. 

3) Changing the relaxation with one having better smoothing properties, of GS type, 
showed large improvements in computer experiments. 

4 ) Preconditioning of the transfers can be used to reduce the amplitudes of the oscillatory 
components. For example smoothing can be introduced before the transfer, say using a sten- 
cil (12 1 )/4 which multiplies the Fourier components with the factor (1 + cos(nk/m))/2. 

5) Higher order interpolations can be used to reduce the Gk,k+m terms. 

6) The algorithm can be optimized using a subroutine which computes the matrix W, for 
different parameters such as A, a, number of relaxation, coarse level relaxation type. An 
optimization search for parameter combinations providing an efficient multilevel cycle can 
be easily performed having such a subroutine. An observation which reduces much the 
analysis of W is that IF has the desired eigenvector (1,0, ....0) T iff its first column is a 
multiple of this vector. An optimization approach can be directed towards the treatment 
of the first column. 

7) An analysis similar with the one performed for the first vector can be performed for 
the first cluster. The elements below the block corresponding to this cluster should be 0 to 
avoid mixing with other clusters. The block corresponding to the cluster may have nonzero 
subdiagonal elements. This suggests that a separation inside the cluster is required. This 
separation can be performed on coarse or on fine levels. 

8) The analysis of W can show what clusters have to be completed and which components 
have to be treated simultaneously. This is important for robustness and efficiency. It is 
simple to observe that the algorithm can be very efficient in converging to a complete clus- 
ter, treated simultaneously, while it will fail if it will treat only one of the components. A 
comparison between a simultaneous and a sequential algorithm can be performed using W. 

9) If a cluster mixes with a second cluster it is an indication to include the second cluster 
in computation and to treat it simultaneously with the first one. 

10 ) A complete cluster, w hich does not mix with another cluster is a good basis for a stable 
subspace technique. 

11 ) Divergence of some components does not imply that the algorithm is not efficient. Sep- 
aration is important. For example, for finding the second cluster one may have to use in 
computations the first cluster too. The first cluster may diverge but if the tw r o clusters get 
well separated from the remaining components, then an algorithm treating both clusters 
simultaneously should be efficient . 
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12 ) Divergence can be treated using different relaxations, e.g., which amplify most speci- 
fied components. 

13 ) The computation of W and the optimization of the algorithm using W can be efficiently 
performed on coarse levels. The computation of W on fine levels is usualy a very expensive 
task but generally not needed, a coarse version of W being sufficient for optimization. 

14 ) The relative separation factors may be important in an FMG algorithm, and not as 
much the asymptotic factors. In an FMG algorithm for computing U \ , the amplitudes B 
come already close to (1, 0, so that only the relative separation of certain components 
may be relevant, usually the ones which mix with U \ . 

15 ) The separation analysis can be useful also for cases when W is not linear. Two such 
cases are when a projection is used in the algorithm, and when W is a composition of 
separation of operators some of which can be analysed. 

16 ) The computation of IT can be performed in different ways. The presented way, i.e., 
using its analytic form, may not be the best one. A way to compute easily IT for different 
changes of the algorithm is required. One way is to compute directly W using the action 
of the algorithm on a basis, e.g., on the columns of U or on the columns of the identity 
matrix. Another way is to compute IT as a composition of simpler mappings. 

17 ) The optimization can be performed by an MG procedure. This is an optimization 
problem with several local minima usually. See [11] for an outline and application of an 
MG optimization approach usable for several local minima. 

18 ) The optimization of the algorithm using W leads to robust algorithms. The asymptotic 
convergence rates can be accurately predicted. Worst cases can be found. 

The behavior predicted by W was accurately reproduced by a program implementing 
the two level algorithm. 
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